Bifurcation and chaos in zero Prandtl number convection 
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We present the detailed bifurcation structure and associated flow patterns near the onset of zero 
Prandtl number Rayleigh Benard convection. We employ both direct numerical simulation and a 
low-dimensional model ensuring qualitative agreement between the two. Various flow patterns origi- 
nate from a stationary square observed at a higher Rayleigh number through a series of bifurcations 
starting from a pitchfork followed by a Hopf and finally a homoclinic bifurcation as the Rayleigh 
number is reduced to the critical value. Global chaos, intermittency, and crises are observed near 
the onset. 

PACS numbers: 47.20.Bp, 47.20.Ky, 47.27.ed, 47.52. +j 



Thermal convection is observed almost everywhere in 
the universe: industrial appliances, liquid metals, atmo- 
sphere, oceans, interiors of planets and stars, galaxies 
etc. An idealized version of convection called Rayleigh 
Benard convection (RBC) has been studied for almost 
a century and it is still an area of intense research 
The two most important parameters characterizing con- 
vection in RBC are the Rayleigh number, describing the 
vigour of buoyancy, and the Prandtl number, being the 
ratio of kinetic viscosity and thermal diffusivity. Solar [2j 
and geological flows [3j are considered to have very low 
Prandtl numbers, as do flows of liquid metals [4|. RBC 
exhibits a wide range of phenomena including instabili- 
ties, patterns, chaos, spatio-temporal chaos, and turbu- 
lence for different ranges of Rayleigh number and Prandtl 
number [f]. Origin of instabilities, chaos, and turbulence 
in convection is one of the major research topics of con- 
vection. 

Direct numerical simulation (DNS), due to its high 
dimensionality, generates realistic but excessively volu- 
minous numerical outputs which obscure the underlying 
dynamics. Lower dimensional projections lead to mod- 
els which, if done improperly, lose the overall physics. 
In this letter, our aim is to unfold and discover the un- 
derlying physics of low Prandtl number flows [f| by ex- 
amining the natural limit of zero Prandtl number (zero- 
P) [1, H S US, EL G3- This offers a dramatic simpli- 
fication without sacrificing significant physics, as well as 
displays a fascinatingly rich dynamic behaviour. In par- 
ticular, since zero-P flows are chaotic immediately upon 
initiation of convection, we adopt a nonstandard strat- 
egy of approaching this system from the post-bifurcation 
direction. Moreover, we attack the problem simultane- 
ously with DNS (to ensure accuracy) as well as a low 
dimensional model (to aid physical interpretation); and 
we stringently refine both the model and DNS until sat- 
isfactory agreement is obtained at all levels of observed 
behaviour. Our results show a diverse variety of both 
new and previously observed flow patterns. These flow 
patterns emerge as a consequence of various bifurcations 
ranging from pitchfork, Hopf and homoclinic bifurcations 
to bifurcations involving double zero eigenvalues. 



Convection in an arbitrary geometry is quite complex, 
so researchers have focused on convective flow between 
two conducting parallel plates called Rayleigh Benard 
convection fl|. The fluid has kinematic viscosity v, ther- 
mal diffusivity k, and coefficient of volume expansion a. 
The top and bottom plates are separated by distance d 
and are maintained at temperatures T\ and Ti respec- 
tively with T\ > T2- Convective flow in RBC is charac- 
terized by the Rayleigh number R = a{T\ — T-i)gd? jvn, 
where g is the acceleration due to gravity, and the 
Prandtl number P = vjn. Various instabilities, pat- 
terns, and chaos, are observed for different ranges of R 
and P 0, H, [HI]. Transition to chaotic states through 
various routes have been observed in convection fLU [la ] . 

In this letter, we focus on zero-P convection. The gov- 
erning zero-P Boussinesq equations 0] are nondimension- 
alized using d as length scale, d 2 jv as time scale, and 
f (Ti — T 2 )/k temperature scale which yields 

d t (W 2 v 3 ) = V 4 v 3 + RV 2 H 6 

-e 3 -Vx[(u-V)v-(vV)u], (1) 

d t u 3 = V 2 w 3 + [(w-Vh-(v-VH, (2) 

V 2 # = -v 3 , (3) 

V-v = 0, (4) 

where v = Va, v 3 ) is the velocity field, 9 is the devia- 
tion in the temperature field from the steady conduction 
profile, u> = V x v is the vorticity field, e 3 is vertically di- 
rected unit vector, and V 2 ^ = d xx + d vy is the horizontal 
Laplacian. We consider perfectly conducting and free- 
slip boundary conditions at the top and bottom plates, 
and periodic boundary conditions along horizontal direc- 
tions @, [HI . In the following discussion we also use the 
reduced Rayleigh number r = R/R c where R c is the crit- 
ical Rayleigh number. 

Straight two-dimensional (2D) rolls that have zero vor- 
ticity are neutrally stable solution of zero-P convection 
at r = 1. However they become unstable for r > 1. 
Busse @, Thual [1] and Kumar et. al. Q showed that 
these 2D rolls saturate through generation of vorticity 
(wavy rolls) for r > 1 both for low Prandtl number and 
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zero-P fluids. Thus vorticity plays a critical role in zero-P 
convection. 

Herring [Io| was first to simulate these equations under 
free-slip boundary conditions. However he observed di- 
vergence of the solutions possibly due to the instabilities 
described above. The first successful simulation of zero-P 
equations with free-slip boundary conditions was done by 
Thual Q . He reported many interesting flow patterns in- 
cluding relaxation oscillation of square patterns (SQOR) 
and stationary square patterns (SQ). Later Knobloch [11| 
explained the stability of the SQ patterns using the am- 
plitude equations. Pal and Kumar [l2| explained the 
mechanism of selection of square patterns using a 15- 
dimensional model. Note that asymmetric squares (re- 
ferred to as 'cross roll' in literature) and other patterns 
have been observed in experiments of low-Prantl number 
convection [H. ITrjl]. 

We performed around 100 DNS runs of zero-P convec- 
tion (Eqs. JT][4])) using a pseudo-spectral code for various 
r values on 64 3 box. The aspect ratio of our simulation 
is 2\/2 : 2y2 : 1. In DNS we observe stationary squares 
(SQ), stationary asymmetric squares (ASQ), oscillatory 
asymmetric squares (OASQ), relaxation oscillations with 
squares (SQOR), and chaos. 

For our bifurcation analysis we construct a low- 
dimensional model using the energetic modes of the 
above-mentioned simulation in the range of r = 1 — 1.4. 
We pick 9 large-scale vertical velocity modes (real Fourier 
amplitudes): W W i, Won, W 20 2, W 02 2, W103, W i3, Wu 2 , 
W\2\, W211, and 4 large-scale vertical vorticity modes 
(real Fourier amplitudes): Zuo, Z112, Z121, Z211. The 
three subscripts are the indices of wavenumber along 
x,y, and z directions. Cumulative energy contained in 
these modes ranges from 85% to 98% of the total en- 
ergy of DNS, and each of these modes has 1% or more of 
the total energy. We derive the model equations by the 
Galerkin projection of Eqs. {T][4]) on the subspace of these 
modes. This results in thirteen coupled first-order ordi- 
nary differential equations for the above variables. The 
low-dimensional model captures all the flow patterns of 
DNS mentioned above. The range of r for these patterns 
for the model and DNS are shown in Table HI and they 
are reasonably close to each other. Interestingly, the sta- 
ble steady values of the modes W101, Won, W112, W121, 
Wm for SQ and ASQ patterns match with correspond- 
ing DNS values within 10%. 

The origin of the above flow patterns can be nicely 
understood using the bifurcation diagram of the low- 
dimensional model. To generate the bifurcation diagram, 
we first evaluate a fixed point numerically using Newton- 
Raphson method for a given r. The branch of fixed 
points is subsequently obtained using a fixed arc-length 
based continuation scheme [l3]- Stability of the fixed 
points is ascertained through an eigenvalue analysis of 
the Jacobian and accordingly the bifurcation points are 
located. New branches of fixed points are born when 
the eigenvalue (s) become zero (pitchfork), and branches 
of periodic solutions appear when the eigenvalue (s) be- 



Flow patterns 


r (Model) 


r (DNS) 


Chaotic 


1 - 1.0045 


1 - 1.0048 


SQOR 


1.0045 - 1.0175 


1.0048 - 1.0708 


OASQ 


1.0175 - 1.0703 


1.0709 - 1.1315 


ASQ 


1.0703 - 1.2201 


1.1316 - 1.2005 


SQ 


1.2201 - 1.4373 


1.2006 - 1.4297 



TABLE I: Range of reduced Rayleigh number r correspond- 
ing to various flow patterns observed in the model and the 
DNS. Here SQ, ASQ, OASQ, and SQOR represent stationary 
squares, stationary asymmetric squares, oscillatory asymmet- 
ric squares, and relaxation oscillation of squares respectively. 



come purely imaginary (Hopf). Subsequent branches are 
generated by calculating and continuing the new steady 
solutions close to the bifurcation points. 

Fixed points are the backbone of bifurcation diagram. 
For r < 1, the origin is the unique stable fixed point corre- 
sponding to the pure conduction state. There is a double 
zero eigenvalue at r = 1 [18], and all the fixed points (13 
in number) arising from r = 1 are unstable for r > 1. 
These fixed points are shown as dotted lines in Fig. Q] 
Four of these branches of fixed points bifurcate from the 
origin; these fixed points satisfy |Wioi| = |Won|. The 
other 8 branches of unstable fixed points emerge from 
nonzero W101 or Won, and they obey |Wioi| 7^ | Won | 
(see Fig. [I}. With an increase of r, these 8 branches 
become stable and merge with the 4 branches that origi- 
nate from the origin. In Fig.[TJ modes W101 and Won are 
presented even though all other modes are also nonzero. 




FIG. 1: Three dimensional view of the bifurcation diagram 
showing the fixed points with solid and dashed curves rep- 
resenting the stable and unstable fixed points respectively. 
Black, blue, and cyan curves represent stationary squares 
(SQ), asymmetric stationary squares (ASQ), and conduction 
state respectively. All the points on the axis (purple lines) 
are 2D roll solutions. 
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After a discussion on fixed points, we focus on the com- 
plete bifurcation diagram shown in Fig. [5] Chaotic solu- 
tions are observed at the onset of convection itself, i.e., 
just above r = 1. A better insight into the origin of the 
various solutions is facilitated by starting the analysis 
at a higher r value and tracking the various bifurcations 
while approaching r = 1. 




r 

FIG. 2: Bifurcation diagram of the model for 0.95 < r < 1.4. 
Fhe stable branches corresponding to stationary squares (SQ) 
and stationary asymmetric squares (ASQ) are represented by 
solid black and solid blue lines respectively. Red, green, and 
brown curves represent the extrema of oscillatory asymmet- 
ric squares (OASQ), relaxation oscillation of squares (SQOR), 
and chaotic solutions respectively. A zoomed view of the bi- 
furcation diagram for the chaotic regime is shown in the inset. 
Branches corresponding to the unstable fixed points are rep- 
resented by dashed lines. Cyan line represents the conduction 
state. 

We start our analysis at r = 1.4 where we observe 
stable symmetric squares (SQ) with |Wioi| = |Woii| 
(black curve in Fig. [T]). In Fig. [2] we represent only the 
Wioi = Won solution. As r is reduced from 1.4, the 
SQ branch of fixed points loses stability via a supercriti- 
cal pitchfork bifurcation at r « 1.2201, after which we 
observe stationary solutions with Wioi 7^ Won (blue 
curves of Figs. Q] and [5]). These solutions correspond 
to asymmetric square patterns (ASQ), either dominant 
along x axis (|Wioi| > |Woii|), or dominant along y axis 
(|Wioi| < |Won|). The SQ solution |W 101 | = |W n| 
continues as a saddle. With a further reduction of r, 
ASQ branches lose stability through a supercritical Hopf 
bifurcation at r w 1.0703 and limit cycles are born. 
These limit cycles are represented by red curves in Fig. [5] 
Physically they represent oscillatory asymmetric square 
patterns (OASQ). Fig. [^a) illustrates the projection of 
two of these stable limit cycles (for r = 1.0494) on the 
Wioi — Won plane. 

The limit cycles grow in size as r is lowered. A homo- 
clinic orbit is formed at r » 1.0175. Afterwards, homo- 
clinic chaos is observed in a narrow window. At lower r 
the attractor becomes regular resulting in a larger limit 
cycle that corresponds to the relaxation oscillations with 
an intermediate square pattern (SQOR). Fig. [3J&) illus- 
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FIG. 3: Phase space projections of limit cycles on the VF101 — 
Won plane for (a) r = 1.0494 and (b) r = 1.0099. The limit 
cycles in figure (a) merge to form a single limit cycle in figure 
(b). Black dots indicate the symmetric square saddle. 



trates the projection of this limit cycle at r = 1.0099. 
The flow pattern in this regime changes from an approx- 
imate pure roll in one direction to a symmetric square, 
and then to an approximate pure roll in the perpendicu- 
lar direction. The SQOR solution is represented by the 
green curve in Fig. [2] 

The flow becomes chaotic as r — ► 1. The chaotic flow 
manifests itself in three different forms: Chi, Ch2, and 
Ch3 as shown in the inset of Fig. [5] as (i), (ii), and (iii) 
respectively. The phase space projection for these three 
solutions are depicted in Fig. 0] for r = 1.0041, 1.0038 
and 1.0030 for the 13-mode model, and for r = 1.0045, 
1.0032 and 1.0023 in the DNS. Their chaotic nature is 
confirmed by the positivity of the largest Lyapunov ex- 
ponent (0.0131, 0.0254 and 0.0389) calculated using the 
13-mode model for r = 1.0041, 1.0038 and 1.0030 respec- 
tively. 

The first chaotic solution, Chi (Figs. H^a) andg^b)) 
results from the broadening of the limit cycle attractor 
with chaotic switchings between the two lobes of the at- 
tractor. This global chaos could probably be attributed 
to homoclinic tangles. The time-series of the solution 
shows intermittency. At r 1.004, these four chaotic at- 
tractors of Chi merge in a 'crisis' to yield a single large 
chaotic attractor Ch2 (Figs.QJc) and0Jd)). This chaotic 
solution persists till r « 1.0035 after which it splits into 
four different chaotic attractors Ch3 in another 'crisis', 
one of which is shown in Figs. [He) and[3Jf). The time- 
series again shows intermittency, and the flow pattern 
switches from an approximately pure roll in one direc- 
tion to an asymmetric square. With a further reduc- 
tion in the Rayleigh number, the size of these chaotic 
attractors decreases and they ultimately merge with one 
of the branches of the unstable ASQ fixed points at 
r = 1 . In Fig. [21 we exhibit the merger of one of these 
chaotic attractors with the unstable ASQ fixed point with 

Wioi -> o. 

In conclusion, we present for the first time a numeri- 
cally obtained, DNS validated, detailed bifurcation dia- 
gram and associated flow structures of zero-P convective 
flow near the onset of convection. The whole spectrum of 
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FIG. 4: The three different chaotic solutions observed near 
r = 1: Chi at r = 1.0041 for the model (a) and at r = 1.0045 
in DNS (b); Ch2 at r = 1.0038 for the model (c) and at 
r = 1.0030 in DNS (d); Ch3 at r = 1.0030 for the model (e) 
and at r = 1.0023 in DNS (f). These solutions belong to (i), 
(ii), and (iii) regimes in the bifurcation diagram (Fig. [2}. 



phenomena observed in DNS near the onset of convection 
is replicated by the low-dimensional model. Hence, the 
bifurcation structure presented here explains the origin 
and dynamics of various patterns near the onset of con- 
vection. Recent analysis of VKS (Von-Karman-Sodium) 
experimental results indicate a strong role of large-scale 
modes for the magnetic field reversal A study of 

large-scale modes as outlined in this letter may provide 
useful insights into the mechanism behind the generation 
and reversal of magnetic field. The dynamics of large- 
scale modes in other hydrodynamic systems like rotating 
turbulence, magneto-convection etc. could also be cap- 
tured by a similar approach. 

A careful analysis of DNS results indicate the existence 
of fringe attractors apart from the main attractors pre- 
sented in this letter. These new attractors are relatively 
insignificant and they exist only in localized regimes of r. 
A detailed investigation of all the attractors will be re- 
ported elsewhere. Bifurcation analysis for r > 1.4 is rea- 
sonably complex as well. Also preliminary results show 
a reasonable amount of similarity between low Prandtl 
number convection and zero-P convection. These issues 
are under investigation. 
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